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We study the thermodynamics of ultracold Bose atoms in optical lattices by numerically diago- 
nalizing the mean-field Hamiltonian of the Bose-Hubbard model. This method well describes the 
behavior of long-range correlations and therefore is valid deep in the superfluid phase. For the 
homogeneous Bose-Hubbard model, we draw the finite-temperature phase diagram and calculate 
the superfluid density at unity filling. We evaluate the finite-temperature effects in a recent ex- 
periment probing number fluctuation [Phys. Rev. Lett. 96, 090401 (2006)], and find that our 
finite-temperature curves give a better fitting to the experimental data, implying non-negligible 
temperature effects in this experiment. 

PACS numbers: 03.75.Lm,03.75.Hh,67.40.-w 



5h 

43 



c 

o 
o 



> 
o 
cn 

(N 
On ■ 

o ■ 
o . 

-I— > ' 

a ■ 

: 

o . 
o ■ 



X 



I. INTRODUCTION 

The ultracold Bose atoms in optical lattices have 
opened a new window to investigate the strongly cor- 
related systems with highly tunable parameters [l( . The 
basic physics of these ultracold atomic systems is cap- 
tured by the Bose-Hubbard model, whose most funda- 
mental feature is the existence of superfluid to Mott- 
insulator (MI) phase transition [1, H| . In a shallow op- 
tical lattice, the ultracold atoms are in superfluid phase, 
which can be well-described by a macroscopic wave func- 
tion with long-range phase coherence Q. In this case, 
phase fluctuation vanishes and on-site number fluctua- 
tion diverges. Whereas, in a deep optical lattice, the 
atoms are in MI phase if the filling factor is fixed to an 
integer, leading to zero on-site number fluctuation and 
divergent phase fluctuation 0, 0] . 

The physics of the MI phase is that, when the repulsive 
interaction between atoms is large enough, the number 
fluctuation would become energetically unfavorable, forc- 
ing the system into a number-squeezed state. This inter- 
action induced MI phase plays an important role in the 
strongly correlated systems, as well as various quantum 
information processing schemes [f|- In the past, a se- 
ries of ultracold-atom experiments have been performed 
to detect this number-squeezed MI phase through the 
observation of increased phase fluctuations [3, M, 7| or 
through an increased time scale for phase diffusion 8 , ]j| . 

Recently, a continuous suppression of on-site number 
fluctuation was directly observed by Gerbier et al. by 
monitoring the suppression of spin-changing collisions 
across the superfluid to Mott-insulator transition [9( . By 
using a far off-resonant microwave field, the spin oscilla- 
tions for doubly occupied sites can be tuned into reso- 
nance, and its amplitude is directly related to the prob- 
ability of finding atom pairs in a lattice site. It was 
shown by Gerbier et al. that, for a small atom num- 
ber, the oscillation amplitude is increasingly suppressed 
when the lattice depth is increased and completely van- 
ishes for sufficiently large lattice depth. They also com- 
pared the experimental results with the prediction of the 



Bose-Hubbard model within a mean-field approximation 
at zero temperature. However, the theoretical curves at 
zero temperature do not fully fit with the experimental 
data, especially in low optical lattices. An important rea- 
son for this deviation is the neglect of finite-temperature 
effects in the calculation. This motivates us to study the 
finite-temperature properties of these ultracold atoms ex- 
pecting a better fitting to the experimental data. 

In the past, various theoretical approaches have been 
used to investigate the zero temperature properties of 
the Bose-Hubbard model such as the mean-field 



proximation 0, fToL EH , strong-coupling expansion 
Gutzwiller projection ansatz [3|, Il3l . Il4j . and quantum 
Monte Carlo simulation [l5(. There are also a few works 
focusing on the finite-temperature properties of the Bose- 
Hubbard model, such as the standard basis operator 
method [HI EH] , slave particle approach [13, EH , coarse- 
graining mean- field approximation [l9l |. analytic inves- 
tigation of fixed point [2(|, and some others [U |22| . 
However, most of these approaches are based on the per- 
turbation theory in terms of small superfluid order pa- 
rameter and therefore are only valid near the superfluid 
to normal-liquid phase transition or in the Mott region. 
Recently, Oosten et al. have used an effective approach, 
which is a finite-temperature extension of the zero tem- 
perature mean field theory proposed by Sheshardi et al. 
[l(|, to study the thermodynamics of the system with 
large filling factor [23[. This mean- field theory correctly 
includes the long-rang correlations fldj . [24[ and thus is 
valid deep in the superfluid phase [25|. In this paper, 
we will use it to calculate the thermodynamic quantities 
from the high superfluid region to the normal-liquid (or 
MI) region. 

This paper is organized as follows. In Sec. [Til we 
will describe the finite-temperature mean-field theory 
for the Bose-Hubbard model. We will draw the finite- 
temperature phase diagram and investigate the system 
with unity filling factor. In Sec. IIII1 we will calculate 
the finite-temperature number fluctuations and compare 
them with the experimental results of Ref. In Sec. 
IIV1 we will give our conclusions. 
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II. FINITE-TEMPERATURE MEAN-FIELD 
THEORY TO BOSE-HUBBARD MODEL 

We consider an ultracold Bose atom gas confined in a 
three-dimensional optical lattice potential. In real exper- 
iments, a slow varying trapping potential is superimposed 
onto the lattice potential. We only pay attention to the 
homogeneous case in this section and will consider the 
trapping potential in the next section. The homogenous 
Bose atom system can be well-described by the following 
Bose-Hubbard Hamiltonian [3[: 

H=-t + -g (!) 

<ij> i i 

Here b\ is the creation operator at site i, nt = b\bi is the 
particle number operator, and (ij) denotes the sum over 
nearest neighbor sites, t and U are the hopping ampli- 
tude and on-site interaction, respectively. In the mean- 
field approximation [10|, the hopping term is decoupled 
as 



b% 



(2) 



where <p = (b\) = is the superfluid order parame- 
ter. The Hamiltonian of the Bose-Hubbard model can 
be written as a sum over single-site terms, H — Hi, 
where 



Hi 



U 



ni(m - I) - firii - zt<p(b] + bi) + zt(f> 2 , (3) 



with z being the number of nearest neighbors. In Ref. 
[icj |. Sheshardi et al. have studied the zero-temperature 
properties of the Bose-Hubbard model by diagonalizing 
the Hamiltonian Hi in the occupation number basis {\n}} 
truncated at a finite value rit- Here, we will extend this 
method to include the temperature effects [23[ and in- 
vestigate the finite-temperature properties of the same 
model. We first obtain the matrix of Hi in the trun- 
cated occupation number basis, which has a symmetric 
tridiagonal form: 



(d{l) e(l) 

e(l) d(2) e(2) 
e(2) d(3) 



V o 







\ 



d{n t -l) e(nt-l) 
e[n t - 1) d{n t ) e(n t ) 
e(n t ) d{n t + 1) / 



e(k) — —Vkzt(j), 



(6) 



and all other matrix elements are zero. We diagonal- 
ize this matrix to obtain the energy spectrum {E^\ and 
eigenstates \ipk), and then evaluate the partition function 
and the free energy, 



k=l 



F 



1 

1 



■ In Z. 



(7) 



For given U , t, p, and T , the superfluid order parameter 
<j> can be determined by minimizing the free energy, i.e., 



dF 



= 0. 



(8) 



U,t,fi,T 



The region with nonzero (f> is identified as the superfluid 
phase while the region with 4> = as the Mott-insulator 
or normal-liquid phase. After determining </>, it is easy to 
calculate other physical quantities such as the superfluid 
density p s and average density p with 



Ps = 



p= (n) = 



Tr(ne-P H ) 



(9) 



,(4) 



where the diagonal elements d(k) and subdiagonal ele- 
ments e(k) are 



d{k) 



U 



(k - l)(Jfc - 2) - p(k - I) + zt(/) 2 



(5) 



We show our main results in Figs. Q] - 2J The finite- 
temperature phase diagrams are plotted in Fig. Q] in 
the p/zt-U/zt plane. The different curves represent the 
phase boundaries between superfluid and normal liquid 
(or MI) at different temperatures. One can see that the 
zero-temperature Mott lobes are enlarged and disappear 
gradually with the increasing temperature. This is the 
so-called MI to normal-liquid crossover. We see that the 
crossover temperature is about 0.5 — §.7zt around the 
first Mott lobe (red curves), while it is about 1.1 — l.izt 
around the second one (blue curves). Note that we set 
fee = 1 throughout this work. The second Mott lobe is 
more stable against the temperature than the first be- 
cause of the larger interaction U. 

In Fig. [21 we show the average density p as a function 
of p/zt for different temperatures and fixed interaction 
U / zt = 12. At zero temperature, there are plateaus on 
the curve with integer filling factor and zero compressibil- 
ity. When increasing the temperature, these plateaus are 
shrunk and eventually disappear. The qualitative behav- 
iors of all these curves are the same as those of previous 
finite-temperature works based on perturbation theory 
[H, E3) O EH- However, our results persist deep into 
the superfluid phase. 

We then turn to investigate the system with an integer 
filling factor, p = 1, at which the superfluid to Mott- 
insulator phase transition may occur. The evolution of 
superfluid density as a function of interaction at different 
temperatures is shown in Fig. [3] At zero temperature, 
p s drops to zero at a critical interaction U c /zt = 5.83, 
in good agreement with the analytical result [ll| . When 
increasing the temperature, the superfluid density as well 
as the critical interaction is reduced. Both of them vanish 
when T / zt is larger than 1.45. We record the critical 
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FIG. 1: (Color online) The finite-temperature phase diagram 
of the Bose-Hubbard model in the fj,/zt-U/zt plane. The 
temperature is increasing when evolving away from the zero- 
temperature Mott lobes. The temperature T/zt of each curve 
is listed in a column on the right side of the figure. The 
red curves with values 0.5, 0.6, and 0.7 on them indicate the 
crossover region for the first Mott lobe, and the blue ones 
with 1.1, 1.2, and 1.3 on them indicate the crossover for the 
second lobe. 




FIG. 2: Average density p as a function of pi/zt at a fixed 
interaction U/zt = 12. The temperature is increasing when 
evolving away from the zero-temperature curve, which is in- 
dicated by the arrows in the figure. We list the temperature 
of each curve in the column on the right side of the figure. 



interaction for every temperature, and plot inversely the 
critical temperature T c /zt as a function of interaction in 
the inset. This is consistent with our previous result in 
Ref. [3. 

In Fig. |H we show the probability P(2) of a site to be 
doubly occupied as a function of interaction U / zt for dif- 
ferent temperatures. At zero temperature, there is an ob- 
vious suppression of P(2) with the increasing interaction, 
i.e., P(2) reduces to zero when entering the MI region. 
However, one can see that the suppression is weakened 
as the temperature is increased. This phenomenon has 
recently been observed by Gerbier et al. by using a spin- 



FIG. 3: Superfluid density p s vs interaction U/zt for different 
temperatures in the case of integer filling p = 1. The temper- 
ature of each curve is increasing from up to down, with their 
values listed in the column on the right side. In the inset, the 
critical temperature T c /zt of the phase transition is plotted 
as a function of interaction U / zt. 
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FIG. 4: Probability P(2) of a site to be doubly occupied vs 
interaction U/zt for different temperatures, with fixing filling 
factor p = 1. The temperature of each curve is decreasing 
from up to down, with their values listed on the right side as 
usual. 



changing collision technique Q. We will discuss this in 
detail in the next section. 

In the final of this section, we evaluate the effect of 
finite truncation to the occupation number basis. As we 
know, the smaller U/zt the larger nt is needed for small 
error |10j. We then plot in Fig. [5] the superfluid densities 
p s of a unity filling system versus the truncation nt for 
a very small interaction U/zt = 0.2. One can see that 
the finite-truncation approximation is quite worse in the 
small n t region, i.e., the magnitude of p s with n t is quite 
different from that with n t + 1. This finite-truncation 
effect is strengthened when increasing the temperature. 
However, when n t is larger than 6, it becomes very small 
even for a high temperature. Typically, the difference of 
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We now specify the values of the parameters used in 
our calculations. The hopping amplitude t and on-site 
interaction U of the Bose-Hubbard model can be deter- 
mined by 



-^V 2 + Vol (r))™(r-r ,),(13) 



t = drw*{v — Yi 



(14) 



where a s is the s-wave scattering length and w(r) is the 
Wannicr function. Here, we use the approximate expres- 
sions I26I 



FIG. 5: Superfluid density p s vs truncation n t for different 
temperatures. The density is p = 1 and the interaction is 
U/zt = 0.2. The temperature of each curve is increasing from 
up to down, with their values listed on the right. 



p s between n t — 10 and n t = 11 is smaller than 0.5%, 
which means the effect of finite truncation is very small. 
Therefore we choose nt = 10 in all of our calculations. 



III. THE FINITE-TEMPERATURE NUMBER 
FLUCTUATIONS: COMPARING TO 
EXPERIMENT 

In this section, we will evaluate the finite-temperature 
number fluctuations of ultracold atoms in optical lat- 
tices and compare them with the experimental data in 
Ref . @ . The ultracold atomic gas is confined in a three- 
dimensional optical lattice potential Vql(j) with 



t 



3 

V OL (r) = Vo5>in 2 (fcr,), 



(10) 



where Vq is the lattice depth, k = 2ir/\ is the wave 
vector with A being the laser wavelength. In addition, 
the atomic gas is also subjected to a trapping potential 
V ex t(r), which can be approximately considered as a har- 
monic one [2611 , 



V ext (r) 



1 



2 2 

muj r , 



with the trapping frequency given by 



8Vn 



(11) 



(12) 



Here m is the atom mass, u m is the frequency of the 
magnetic trap where the condensate is initially formed, 
and w is the waist of the lattice beams 0, |2f|. If this 
trapping potential is very shallow and varies slowly across 
the atomic cloud, we can treat it within the local density 
approximation (for more details, see Ref. (2r|). 



0.08 



— = 1.43 i ) exp (-2.O7V%/^0> (15) 



V_ 



= 5.97 



(f) 



0.88 



(16) 



where E r — h 2 /2m\ 2 is the single-photon recoil en- 
ergy. The triplet scattering length for the 87 Rb atom 
is a s = 5.61 nm [27|, the laser wavelength A = 842 nm, 
the lattice spacing d = A/2 = 421 nm, the magnetic trap- 
ping frequency uj m = 2ir x 16 Hz, and the waist w = 136 
fim. All these values coincide exactly with those in the 
real experiment 0] . 

We show in Fig. [6] the density distribution of ultra- 
cold atomic gas with total number N — 1.0 x 10 5 in two 
typical optical lattice depths. We can see that, at zero- 
temperature and in the low optical lattice [Fig. [Ha)], the 
density decreases smoothly when leaving the center of the 
trap. However, a Mott-shell structure forms in the deep 
lattice [Fig. E(b)]. When increasing the temperature, 
the atomic cloud is expanding and the Mott-plateaus in 
the deep lattice are melting gradually. Most recently, 
the formation of Mott-shell structure is directly observed 
[28l [20( 1 . The calculation of the finite-temperature effects 
in these experiments is ongoing. 

We now turn to the number fluctuation experiment 
performed by Gerbier et at [9j. In explaining their ex- 
perimental results, the authors used a zero-temperature 
mean-field theory, which gives the correct trend. How- 
ever, they also found some discrepancies, especially in 
the low lattice depths [9j. We find that these discrepan- 
cies can be diminished by taking the finite-temperature 
effect into account. We show the probability P(2) of find- 
ing atoms in the doubly occupied sites, averaged over the 
whole atomic cloud, in Fig. [7] as a function of total atom 
number N for different lattice depths and temperatures. 
The zero temperature curve in each panel coincides with 
the theoretical estimate in Ref. 9j. One can see that 
their deviation from the experimental data, which are de- 
noted by circles in the figure, is quite obvious. In panel 
(a) where the atoms are in superfluid phase, P{2) always 
decreases with the increasing temperature. When com- 
paring with the experimental data, we can see that for 
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FIG. 6: Density distribution of ultracold atom gas with a total 
atom number N = 1.0 x 10 J in different optical lattice depths 
Vo = SE r (a) and Vo = 18E r (b). The zero-temperature 
curves are indicated by arrows in the figure. The temperature 
of each curve is increasing when evolving away from the zero- 
temperature one, with their values listed in the columns on 
the right side. 



small atom number the discrepancy between the zero- 
temperature curve and the experimental data can be re- 
duced by taking the temperature into account. In panel 
(b), we observe that the finite-temperature curve with a 
proper temperature, e.g., the red curve with T/zt = 0.7, 
agrees well with the experimental data in the whole range 
of N . In panel (c) where the atoms are near the super- 
fluid to Mott-insulator transition, the influence of the 
temperature is a little complex. In the small N region, 
P(2) first decreases and then increases with the increas- 
ing temperature. The decrease at the beginning is due to 
the suppression of superfluid density by the temperature, 
and the following increase is caused by the temperature 
exciting more atoms to occupy the eigenstate |2). In the 
large N region, however, the behavior of P(2) is inverse: 
it first increases and then decreases with the increasing 
temperature. Again, we find that the curve with high 
temperature (e.g., the red curve with T/zt = 1.5) fits 
fairly well with the experimental data. In panel (d) where 
the atoms are deep in the Mott region, the amplitude of 
P(2) is improved in the small N region, but is suppressed 
in the large N region, as the temperature increases. This 
behavior is consistent with the experimental data too. 
The investigation of these finite-temperature effects has 
its own interest, for it provides a possible way to esti- 



mate the temperature of the atomic gas [28j. In our 
case, l.Ozt yields several ten nano-Kelvin and the tem- 
perature regime where the experimental data located in 
is 10 1 ~ 10 2 nK, therefore our estimate to temperature 
is quite reasonable. 

Before closing this section, we observe that, in the large 
N region, there are still some inconsistencies between our 
results and the experimental data. Within the deep su- 
perfluid regime [Fig. E[a)], P(2) is always smaller than 
the experimental data in the large N region. A possi- 
ble reason for this discrepancy is the large probability 
for triplet occupation, which is due to the high density 
in the center of the trap (p ~ 3) and the large num- 
ber fluctuation of the system in so low lattice potential. 
In experiment, the presence of large triplet occupation 
may change the spin resonance condition, resulting in a 
considerable contribution of spin oscillations on triply oc- 
cupied sites to the observed oscillation amplitude [9l[30|. 
This makes the experimental measurements larger than 
the calculated P(2) (see also Fig. 4 of Ref. @). Within 
the deep MI regime [Fig. E[d)], our finite-temperature 
consideration does give a correct trend in the large N re- 
gion. However, the reason for the further suppression of 
oscillation amplitude is not uncovered yet, which cannot 
be simply attributed to the temperature effects. 

IV. CONCLUSIONS 

We studied the finite-temperature properties of the 
Bose-Hubbard model by numerically diagonalizing the 
mean-field Hamiltonian on a truncated occupation num- 
ber basis. This method can give the thermodynamics 
deep in the superfluid phase. For the homogenous Bosc- 
Hubbard model, we calculated the finite-temperature 
phase diagram and plotted the density versus chemical 
potential curve. We demonstrated the evolution of zero- 
temperature Mott lobes and Mott plateaus when the 
temperature is increasing. We then specially investigated 
the system with unity filling factor. The superfluid den- 
sity and the probability P(2) of finding a site to be dou- 
bly occupied were plotted as a function of interaction. 
We found that the superfluid density is reduced and the 
suppression of P(2) is weakened as the temperature in- 
creases. In the second part of this paper, we evaluated the 
finite-temperature effects in a recent experiment probing 
the number fluctuation Q. We showed that the calcu- 
lated finite-temperature curves fit the experimental data 
better than the zero-temperature ones. This implies that 
the finite-temperature effects in this experiment are quite 
large. 
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FIG. 7: (Color online) The probability P(2) of finding atoms 
in the doubly occupied sites, averaged over the whole atomic 
cloud, vs the total atom number TV for different lattice depths 
and temperatures. The circles are the experimental results 
reproduced from Ref. [jj. In panels (a) and (b), the temper- 
ature of each curve is increasing from up to down. In panel 
(c), the blue dash-dot curve corresponds to T/zt = 0.0. For 
other curves, the temperature is decreasing from up to down 
in the small TV region, but is increasing from up to down in 
the large TV region. Panel (d) has a similar temperature trend. 
From up to down, the temperature is decreasing in the small 
TV region but is increasing in the large TV region. As usual, 
the value of temperature is listed on the right side of each 
panel. The temperatures of the red curves in panels (a), (b), 
(c), and (d) are T/zt =0.7, 0.7, 1.5, and 4.0, respectively. 
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